load "get_spectrum_normal_significant.ncl"
begin
   idir = (/"/g/data/p66/ars599/obs/month/" ,\
            "/g/data/p66/ars599/work_atm_temp/" /)
   iName = (/"HadISST_sst.1970-2019.n96.nc", \
             "ts_Amon_bx944_piControl_r1i1p1_0001-0100.nc", \
             "ts_Amon_bx944a_piControl_r1i1p1_0001-0100.nc", \
             "ts_Amon_bx944b_piControl_r1i1p1_0001-0100.nc", \
             "ts_Amon_bx944c_piControl_r1i1p1_0001-0100.nc", \
             "ts_Amon_bx944d_piControl_r1i1p1_0001-0100.nc" /)
   im1file =  addfile(idir(1)+iName(1),"r")
   im2file =  addfile(idir(1)+iName(2),"r")
   im3file =  addfile(idir(1)+iName(3),"r")
   im4file =  addfile(idir(1)+iName(4),"r")
   im5file =  addfile(idir(1)+iName(5),"r")
   iofile =  addfile(idir(0)+iName(0),"r")

; --- sst bias ---;
   yrStrt = 1995
   yrLast = 2014
   TIME   = iofile->time
   YYYY   = cd_calendar(TIME,-1)/100                 ; entire file
   iYYYY  = ind(YYYY.ge.yrStrt .and. YYYY.le.yrLast)

   obs_all_sst = iofile->sst(iYYYY,:,:)
   mod_all_sst = im1file->ts(iYYYY,:,:)
   mod2_all_sst = im2file->ts(iYYYY,:,:)
   mod3_all_sst = im3file->ts(iYYYY,:,:)
   mod4_all_sst = im4file->ts(iYYYY,:,:)
   mod5_all_sst = im5file->ts(iYYYY,:,:)

   mod_all_sst = mod_all_sst - 273.15
   mod_all_sst = where(mod_all_sst.le.0, 0, mod_all_sst)

   osst = dim_avg_n_Wrap(obs_all_sst, 0)
   dsst = dim_avg_n_Wrap(mod_all_sst, 0)
   dsst = dsst - osst  ; model - obs sst bias
delete([/TIME,YYYY,iYYYY,  dsst@long_name, dsst@standard_name, dsst@units/])
dsst@long_name = "Model minus Obs SST bias"
printVarSummary(dsst)
; --- indices --- ;
   obs_sst = iofile->sst(:,{-5:5},{190:240})
   mod1_sst = im1file->ts(:,{-5:5},{190:240})
   mod2_sst = im2file->ts(:,{-5:5},{190:240})
   mod3_sst = im3file->ts(:,{-5:5},{190:240})
   mod4_sst = im4file->ts(:,{-5:5},{190:240})
   mod5_sst = im5file->ts(:,{-5:5},{190:240})

   obs_ao_sst = iofile->sst(:,{-20:20},{300:360})
   mod_ao_sst = im1file->ts(:,{-20:20},{300:360})
   obs_ao = dim_avg_n_Wrap( obs_ao_sst ,(/1,2/))
   mod_ao = dim_avg_n_Wrap( mod_ao_sst ,(/1,2/)) - 273.15
   obs_ao = runave_n_Wrap (obs_ao, 12, 1, 0)
   mod_ao = runave_n_Wrap (mod_ao, 12, 1, 0)

copy_VarMeta(mod_ao_sst(:,0,0), mod_ao)
printVarSummary(mod_ao)
obs_ao@long_name = "observation ao index"
mod_ao@long_name = "model ao index"

   obs_n34 = dim_avg_n_Wrap(rmMonAnnCycTLL( obs_sst ),(/1,2/))
   mod1_n34 = dim_avg_n_Wrap(rmMonAnnCycTLL( mod1_sst ),(/1,2/))
   mod2_n34 = dim_avg_n_Wrap(rmMonAnnCycTLL( mod2_sst ),(/1,2/))
   mod3_n34 = dim_avg_n_Wrap(rmMonAnnCycTLL( mod3_sst ),(/1,2/))
   mod4_n34 = dim_avg_n_Wrap(rmMonAnnCycTLL( mod4_sst ),(/1,2/))
   mod5_n34 = dim_avg_n_Wrap(rmMonAnnCycTLL( mod5_sst ),(/1,2/))

   obs_n34 = dtrend(obs_n34,False)
   mod1_n34 = dtrend(mod1_n34,False)
   mod2_n34 = dtrend(mod2_n34,False)
   mod3_n34 = dtrend(mod3_n34,False)
   mod4_n34 = dtrend(mod4_n34,False)
   mod5_n34 = dtrend(mod5_n34,False)

   obs_std = dim_stddev( obs_n34 )
   obs_n34 = obs_n34 / obs_std
   mod_n34 = mod1_n34
   mod_n34 = (mod1_n34 + mod2_n34 + mod3_n34 + mod4_n34 + mod5_n34)/5
   mod_std = dim_stddev( mod_n34 )
   mod_n34 = mod_n34 / mod_std
delete([/obs_n34@standard_name, mod_n34@standard_name/])
obs_n34@long_name = "observation N34 index"
mod_n34@long_name = "model N34 index"
printVarSummary(obs_n34)
printVarSummary(mod_n34)
;--- output variables ---
  diro    = "./"
  pltName = "fig01"
  outFile = diro+pltName+".nc"
  system("\rm "+outFile)
  fout         = addfile( outFile,"c")
  fout->obs_ao = obs_ao
  fout->mod_ao = mod_ao
  fout->obs_n34 = obs_n34
  fout->mod_n34 = mod_n34
  fout->dsst = dsst

;--- time axial ---;
  time  = obs_n34&time
  tt    = cd_calendar(time, 4)
  wks_time = tt

;*** ploting spectrum ********************
  optPlt          = True
  optPlt@trYMinF = 0.
  optPlt@trXMinF = 0.1
  optPlt@vpHeightF= 0.1          ; change aspect ratio of plot
  optPlt@vpWidthF = 0.9
  optPlt@prdGap   = 1./12
  optPlt@trXLog   = False
  optPlt@frqOrPrd = False
  optPlt@vpHeightF= 0.4                    ; change aspect ratio of plot
  optPlt@vpWidthF = 0.8
;  optPlt@tiMainString = "Power Spectrum of Nino3.4"
  optPlt@units = "Years"
  pltType  = "png"
  pltName  = "fig01_n34_spec_pattern"

  res                    = True              ; plot mods desired
  res@gsnDraw            = False             ; don't draw yet
  res@gsnFrame           = False             ; don't advance frame yet
  res@vpHeightF 	 = 0.4               ; change aspect ratio of plot
  res@vpWidthF 	         = 0.8
  res@trYMaxF            = 28.0
  res@trYMinF            = 25.6
  res@trXMinF	         = 1970              ; set x-axis minimum
  res@xyMonoLineColor    = False             ; want colored lines
  res@xyLineColors       = (/"Black","Red","Blue"/) ; colors chosen
  res@xyLineThicknesses	 = (/4.,4.,4./)      ; line thicknesses
  res@xyDashPatterns	 = (/0.,0.,0./)      ; make all lines solid

  res@tmXMajorGrid                 = True   ; Turn on vertical lines
  res@tmXMajorGridThicknessF       = 0.5
  res@tmXMajorGridThicknessF       = 0.5
  res@tmXMajorGridLineDashPattern  = 2
  res@tmYMajorGrid                 = True   ; Turn on horizontal lines
  res@tmYMajorGridThicknessF       = 0.5
  res@tmYMajorGridLineDashPattern  = 2

  cnres                       = True     ; plot mods desired
  cnres@gsnDraw               = False    ; don't draw yet
  cnres@gsnFrame              = False    ; don't advance frame yet
  cnres@vpHeightF             = 0.45     ; change aspect ratio of plot
  cnres@vpWidthF              = 0.8 
  cnres@cnFillOn              = True     ; turn on color fill
  cnres@cnLinesOn             = False    ; turn of contour lines
  cnres@cnLevelSpacingF       = 0.5      ; contour spacing
  cnres@cnFillPalette         = "BlueWhiteOrangeRed" ; BlueDarkRed18"
;  cnres@lbOrientation         = "Vertical"
  cnres@mpCenterLonF          = 240.         ; defailt is 0 [GM]
  cnres@gsnAddCyclic          = False    ; data already has cyclic point
				       ; this must also be set for any zoom
  cnres@cnLevelSelectionMode    =  "ManualLevels"
  cnres@cnMinLevelValF          =  -6
  cnres@cnMaxLevelValF          =   6
  cnres@cnLevelSpacingF         =   1
; --- add legend --- ;
  gres = True
  gres@YPosPercent = 90.    ; expressed as %, 0->100, sets position of top border of legend
                            ;  when gres@Position is set to its default setting of "Top" (Default = 95.)
  gres@XPosPercent = 60     ; expressed as %, 0->100, sets position of left border of legend(Default = 5.)
  lineres = True
  lineres@lgLineColors = (/"black"/) ; , "red","orange","dodgerblue"/) ; line colors
  lineres@lgLineThicknesses = 4                          ; line thicknesses
  lineres@LineLengthPercent = 9.                         ; expressed as %, 0->100, length of line
  textres = True
; --- plot --- ;
    wks  = gsn_open_wks(pltType,pltName)       ; specifies a ps plot
    plot = new(4,graphic)
    optPlt@trYMaxF = 40.
    optPlt@trXMaxF = 12.
  optPlt@gsnCenterString = "Simulated Nino3.4 Power Spectrum"
    plot(0) = plot_spectrum_normal(wks,mod_n34,optPlt)              ; raw data
  textres@lgLabels = (/"ACCESS-CM2"/)  ; legend labels (required)
    plot(0) = simple_legend(wks,plot(0),gres,lineres,textres)
  optPlt@gsnCenterString = "Observed Nino3.4 Power Spectrum"
    plot(1) = plot_spectrum_normal(wks,obs_n34,optPlt)              ; raw data
  textres@lgLabels = (/"HadISST"/)  ; legend labels (required)
    plot(1) = simple_legend(wks,plot(1),gres,lineres,textres)

   cnres@tiMainString   = "SST bias"
    plot(2) = gsn_csm_contour_map(wks,dsst, cnres)
   res@gsnCenterString   = "SST over the tropical Atlantic Ocean" ; 20~S~o~N~S-20~S~o~N~N"
   res@tiYAxisString = "~S~o~N~C"                  ; yaxis
    plot(3) = gsn_csm_xy (wks,wks_time,obs_ao,res)  ; Create filled XY plot.
   res@xyLineColors       = (/"Red","Blue","Black"/) ; colors chosen
    plot_a  = gsn_csm_xy (wks,wks_time(0:539),mod_ao(0:539),res)  ; Create filled XY plot.
overlay(plot(3),plot_a)
  textres@lgLabels = (/"HadISST"/)  ; legend labels (required)
    plot(3) = simple_legend(wks,plot(3),gres,lineres,textres)
  gres@YPosPercent = 90.    
delete([/lineres@lgLineColors,textres@lgLabels/])
  lineres@lgLineColors = (/"black", "red"/)
  textres@lgLabels = (/"HadISST", "ACCESS-CM2"/)  ; legend labels (required)
    plot(3) = simple_legend(wks,plot(3),gres,lineres,textres)
; --- panel plot --- ;
    resP = True
;    resP@gsnPanelMainString        = "Power Spectrum of WWB" ; new resource added in NCL V6.4.0
    resP@gsnPanelFigureStringsJust = "TopLeft"                      ; new resource added in NCL V6.4.0
    resP@gsnPanelFigureStrings     = (/"a)","b)","c)","d)"/) ; add strings to panel
    resP@gsnPanelRowSpec = False                   ; tell panel whatorder to plot
    gsn_panel(wks,plot,(/2,2/),resP)
delete(wks)   ; Make sure PS file is closed

;***********************
;---Convert to animated PNG
;***********************
;print("Converting .......")
;    cmd = "convert -geometry 2048 -density 432x432 "+pltName+".eps " +pltName+".png"
;    print(cmd)
;    system(cmd)

end
